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Abstract 

We  discuss  the  application  of  the  reduced  basis  method  for  the  simulation  and  con¬ 
trol  of  unsteady  viscous  flows  governed  by  the  incompressible  Navier-Stokes  equations. 
We  describe  how  to  use  this  method  in  terms  of  the  construction  of  a  lower-order 
compensator  design.  Our  approach  includes  a  construction  method  of  the  optimal 
state  feedback  law  for  finite  dimensional  nonlinear  regulator  problems.  The  method  is 
applied  to  construct  a  feedback  law  for  the  reduced  order  control  model  of  the  Navier- 
Stokes  equations,  and  then  we  apply  our  feedback  law  to  the  original  control  system. 
Our  method  is  demonstrated  on  a  control  problem  formulated  in  a  channel  flow  using 
a  boundary  velocity  control.  We  also  show  how  these  ideas  can  be  extended  to  control 
problems  governed  by  partial  differential  equations.  Numerical  results  are  reported  for 
the  open  and  closed  loop  controls,  and  a  compensator  design  is  proposed  to  complete 
the  closed  loop  dynamics. 

Keywords:  flow  control,  compensator  design,  Navier-Stokes  equations,  reduced  basis 
method,  finite  element  method. 


1  Introduction 

In  this  paper  we  discuss  the  application  of  a  reduced  order  model  for  the  simulation 
and  control  of  unsteady  viscous  flows.  The  reduced  order  model  is  developed  using 
the  reduced  basis  method  with  Lagrange  interpolation.  We  will  use  the  reduced  ba¬ 
sis  method  for  the  construction  of  a  lower  order  compensator  design.  Our  discussion 
will  focus  on  the  active  control  of  the  incompressible  Navier-Stokes  equations,  but  the 
procedure  can  be  applied  to  other  control  problems  governed  by  partial  differential 
equations. 

Active  control  of  partial  differential  equations  in  real  time  requires  extensive  com¬ 
putational  time  and  resources.  Using  standard  methods  such  as  finite  differences  or 
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finite  elements  often  result  in  problems  that  are  too  large  to  solve  in  a  reasonable  time. 
We  develop  a  reduced  order  model  which  overcomes  this  difficulty.  The  reduced  ba¬ 
sis  method  uses  basis  functions  which  are  generated  from  the  equations  being  solved. 
Since  the  basis  functions  are  problem  dependent,  the  dynamics  of  the  original  system 
can  effectively  be  captured  while  using  a  small  number  of  basis  functions.  This  results 
in  a  dense,  low  order  system  which  can  be  solved  efficiently.  This  is  in  contrast  to 
finite  element  methods,  which  use  piecewise  polynomials  as  basis  functions,  resulting 
in  a  large,  sparse  system  of  equations. 

In  §2  we  will  discuss  the  simulation  of  the  incompressible  Navier-Stokes  equations 
using  the  finite  element  method.  In  §3  we  describe  the  reduced  basis  method.  This 
includes  the  selection  of  the  reduced  basis  subspace  and  the  construction  of  the  re¬ 
duced  order  model  for  the  time  dependent  Navier-Stokes  equations.  We  will  compare 
two  techniques  for  the  selection  of  basis  functions,  and  compare  the  reduced  order 
solutions  with  the  full  solution  computed  in  §2.  In  §4  we  discuss  the  reduced  order 
optimal  control  problem  and  formulate  the  solution  method  for  the  nonlinear  regulator 
problem.  In  §5  we  describe  the  construction  of  the  feedback  synthesis  based  on  the 
reduced  order  model,  and  we  also  develop  a  compensator  design  based  on  the  reduced 
order  control  system  and  the  constructed  feedback  synthesis. 


2  Simulation  of  the  Navier-Stokes  Equations 

The  Navier-Stokes  equations  which  govern  viscous  incompressible  flow  are 


— —  Au  +  u  •  Vu  +  Vp  = 

lie 

f  in  fi  X  (0, T] 

(2.1) 

Vu  = 

0  in  0  X  (0,  T] 

(2.2) 

u  = 

U5  on  T  x  [0,  T] 

(2.3) 

u(x,0)  = 

uq(x)  in  0 

(2.4) 

where  u(x,f)  is  the  velocity,  p(x,  t)  is  the  pressure,  f(x,  t)  is  the  body  force,  Re  is  the 
Reynolds  number,  uj  is  the  boundary  velocity,  uo(x)  is  the  initial  velocity,  and  0  is  a 
bounded  region  in  ]R  2  with  boundary  T. 

We  will  cast  (2.1  )-(2.4)  in  a  variational  form  and  use  a  mixed  finite  element  approx¬ 
imation.  First,  we  will  introduce  some  standard  notations.  We  denote  by  i2(fi)  the 
collection  of  square  integrable  functions  defined  on  0,  and  we  define  the  inner  product 
on  L2(  S4)  by 

/  4>V<l-r 

Jn 

Let 

Hl{Q)  =  {u  e  L2(Sl)  :  —  €  i2(fl),  i  =  1,  2} 

H](Sl)  =  {u  e  H^Sl)  :  u\m  =  0} 

The  vector-valued  counterparts  of  these  spaces  are  denoted  by  bold  face  symbols,  for 
example  =  II Mil)  X  //  '($>). 

The  variational  form  of  (2.1)- (2 .4)  is  obtained  by  taking  the  X2-inner  product  of 
Equation  (2.1)  with  any  function  v  €  Hq(I1),  and  the  /Winner  product  of  Equation 
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(2.2)  with  any  function  q  €  i2(fi).  Applying  the  divergence  theorem  we  arrive  at  the 
following  variational  form: 

(ut, v)  +  -^a(u,v)  +  b(  u,  u,  v)  -  c(v,p)  =  (f,  v)  for  all  v  <E  Hq(O)  (2.5) 

c(u,  g)  =  0  for  all  q  E  L2(£l)  (2.6) 

where  (■,•)  denotes  the  i2-inner  product  and  we  define  the  standard  bilinear  and 
trilinear  forms  as  follows: 

a(u,v)  =  /  VuVvdx  for  all  u,  v  6  H1  ( 0 ) 

JQ 

c(u,  q)  =  f  (V  ■  u )qdx  for  all  u  E  H^ff),  q  E  X2(0) 

Jn 

6(u,v,w)  =  -  /  ((u  •  Vv}>  w  —  (u  •  Vw)  •  v) dx  for  all  u,  v,  w  €  H1(fi) 

2  Jn 


where  we  use  the  skew-symmetric  form  of  the  trilinear  form  introduced  in  [1]  so  that 
the  discretized  equations  preserve  the  skew  symmetry  property  of  the  Navier-Stokes 
equations.  We  also  have  the  following  properties  of  these  linear  forms.  All  of  the  forms 
a[  -i  ■),  b(  •,  •.  ■),  and  c(  •,  ■)  are  continuous.  The  bilinear  form  a(  ■.  ■)  satisfies  the  following 
coercive  property: 

a(  u,  u)  >  C'i||u||2  for  all  u  €  Hq(Q). 

The  bilinear  form  c(-,  ■)  satisfies  the  following  inf-sup  condition,  as  explained  in  [2]: 
there  exists  7  >  0  such  that 


inf  sup 
geL2(n)  vgHj(fl) 


(  c(v.g)  A 

V||v||i|k||0J 


>  7- 


(2.7) 


We  will  use  the  mixed  finite  element  method  to  approximate  Equations  (2.5)-(2.6). 
That  is,  we  seek  solutions  E  Vh  and  ph  E  Sh  such  that 

(uf.v'M  +  -^-a(uh,yh)  +  b(uh,  uh,vh)  -  c(vh,ph )  =  (f ,  vh)  for  all  vh  €  (2.8) 

c(u/l,g/l)  =  0  for  all  qh  €  Sh  (2.9) 

where  Vq  C  Hj(fi)  and  Sh  C  L2(Q)  are  approximating  finite  dimensional  subspaces. 
We  use  piecewise  quadratic  and  piecewise  linear  basis  functions  for  the  velocity  and 
pressure,  respectively.  This  selection,  referred  to  as  the  Taylor- Hood  element  pair,  is 
known  to  satisfy  the  condition  (2.7). 

The  resulting  approximating  system  can  be  written  as  the  following  system  of  equa¬ 
tions: 

Qu  +  Hu  +  B(u,  u)  +  =  f  (2.10 ) 

(Bl  B2)  (H)  =  0  (2-U) 

where  Q  is  the  mass  matrix,  H  is  the  stiffness  matrix,  {B\  B2)  are  the  continuity 
matrices,  and  B(u ,  u)  is  the  quadratic  form  that  approximates  the  convective  term.  To 
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solve  this  system,  we  will  project  Equation  (2.10)  onto  the  divergence  free  subspace, 
as  described  in  [3],  by  using  the  following  projection  preconditioner 

P=  PL2(~  A)-1  PL2, 


where 

/  (Jlj  (BiB?  +  B2Bt2  )%!U  lh) 

is  the  X2-projection  onto  the  divergence  free  space.  Applying  the  preconditioner  to 
(2.10  )-(2.11 )  we  have  the  following  nonlinear  system  of  equations 

PQii  +  PHu  +  PB(u,  u)  =  Pf.  (2.12) 


We  used  the  Implicit  Midpoint  rule  to  integrate  Equation  (2.12)  in  time,  ie 


PQ- 


uk+k  -  -i 
At/2 


+  PHuk+1i 


+  PB(uk+12. 


uk+h  =  pj, 


(2.13) 


uk+1  =  2 uk+12  -  uk.  (2.14) 

We  solved  (2.13)-(2.14)  using  Newton’s  method;  the  Newton  step  here  is  carried  out 
using  the  iterative  method  GMRES.  In  Figure  1,  we  see  the  solution  to  these  equations 
on  the  forward  facing  channel  flow.  The  height  of  the  inflow  boundary  is  3,  and  that 
of  the  outflow  boundary  is  2.  For  boundary  conditions,  we  assume  a  parabolic  velocity 
profile  at  the  inflow  and  outflow,  and  zero  boundary  conditions  on  all  other  walls.  The 
inflow  condition  is  u (y)  =  |r/(3  —  y)  and  the  outflow  condition  is  u(j/)  =  §(3  —  y)(y  —  1). 
The  equations  (2.13)-(2.14)  were  solved  using  At  =  0.1  to  a  final  time  of  T  =  40.  Since 


Figure  1:  Time  Dependent  Navier-Stokes  Equations  at  T  =  40 

the  initial  condition  must  also  satisfy  the  incompressibility  condition  (2.9),  we  used 
the  solution  to  the  stationary  Stokes  equations  as  the  initial  condition.  The  Stokes 
equations  are  the  linear  equations  of  the  form 


ut  -  Au  +  Vp  =  0, 
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V  •  u  =  0. 


Typically,  a  large  number  of  unknowns  is  required  to  capture  the  complex  flow 
patterns  at  high  Reynolds  numbers.  For  our  computations,  over  5000  unknowns  were 
used.  Using  an  Ultra  Sun  Sparc  Station,  the  above  computation  required  approximately 
20  hours  of  computational  time.  We  would  like  to  develop  a  method  to  reduce  the  order 
of  the  system,  thus  requiring  less  computational  time  to  solve,  while  still  capturing  the 
dynamics  of  the  original  system.  We  will  describe  this  reduced  order  method  in  the 
next  section. 


3  The  Reduced  Basis  Method 

In  this  section  we  will  present  a  reduced  order  method  based  on  the  reduced  basis 
method.  We  would  like  this  reduced  order  method  to  have  the  following  characteristics: 

1.  Require  a  small  number  of  basis  functions. 

2.  Be  inexpensive  to  solve. 

3.  Capture  the  dynamics  of  the  original  system. 

The  reduced  basis  method  was  first  proposed  in  [4,  5]  for  structural  analysis,  and  it 
has  been  used  for  structural  problems  in  [6,  7],  and  for  incompressible  flow  problems 
in  [8,  9].  We  will  discuss  the  choice  of  reduced  basis  subspace  and  the  construction 
of  the  reduced  order  equation.  We  will  then  develop  a  reduced  order  model  for  the 
incompressible  Navier-Stokes  equations,  and  show  that  it  can  effectively  be  used  to 
capture  unsteady  viscous  flows. 

Consider  the  equation 


£ (u.  A)  =  0  for  ii£l,  A  G  A 


(3.1) 


where  A  represents  a  physical  parameter,  such  as  Reynolds  number,  and  u  is  a  solu¬ 
tion  in  a  Hilbert  space  X .  We  will  construct  the  reduced  basis  subspace  XR  C  X  by 
interpolation  of  the  solution  function  u(  A)  of  Equation  (3.1)  using  one  of  the  following 
three  choices. 

1.  The  Taylor  Subspace  In  this  case  the  reduced  basis  functions  are  linear  combi¬ 
nations  of  Taylor  basis  functions  generated  by  computing  the  Taylor  expansion  of  u{  A) 
at  some  reference  value  of  A,  say  A*.  Let  u0  be  the  solution  at  the  parameter  A*,  and 
assume  it  has  M  derivatives.  The  reduced  basis  subspace  is 


XR 


span{uj  :  uj 


dju 

~dX)  A=A*’ 


j  =  0,...,M} 


The  equations  for  can  be  obtained  by  successive  differentiation  of  Equation  (3.1) 
with  respect  to  the  parameter  A.  For  example,  u\  =  satisfies  the  linear  equation 


£u(u0,  A>i  =  -£a(«o,A*). 


Note  that  each  subsequent  uj .  j  =  I , ....  /I  /  can  be  computed  by  solving  a  linear  sys¬ 
tem  with  the  same  linear  operator  Eu(  u0,  A*).  This  choice  has  been  proven  effective 
for  steady  state  fluid  flow  calculations,  see  [8].  However,  one  disadvantage  is  that  if 
the  parameter  A  moves  away  from  the  reference  parameter  A*,  then  the  basis  functions 
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may  have  to  be  updated. 

2.  The  Lagrange  Subspace  In  this  case  the  reduced  basis  functions  are  linear  com¬ 
binations  of  basis  functions  generated  by  solving  the  nonlinear  system  (3.1)  at  various 
parameter  values  A The  reduced  basis  subspace  is 

XR  =  span{uj  :  Uj  =  -u(  AJ),  j  —  1, . . .  ,M } 


An  advantage  with  this  choice  is  that  the  subspace  can  be  updated  by  adding  more 
basis  functions,  rather  than  updating  the  whole  subspace. 

3.  The  Hermite  Subspace  In  this  case  the  reduced  basis  functions  are  computed  by 
solving  the  nonlinear  system  (3.1)  at  various  parameter  values  AJ,  and  computing  the 
first  derivatives  at  various  parameter  values.  The  reduced  basis  subspace  is 


du , 


X  =  span{uj  =  tt(Aj)  and  —  |A_..A  .  j  =  1, . . . ,  M} 


We  can  think  of  this  choice  as  a  hybrid  of  the  Lagrange  and  Taylor  approaches. 

Let  be  a  basis  for  the  reduced  basis  subspace  XR  C  X ,  and  let  m  =  dim(XR). 
If  we  let  the  reduced  order  solution  be  uR  =  then  we  obtain  the  reduced  or¬ 

der  problem  by  projecting  Equation  (3.1)  onto  the  reduced  basis  space  by  the  Galerkin 
approximation 

Er(ur,  X)i  =  ( E{ur,  A),  dEix-xX--  «  =  !,■■•,  m. 

For  the  evolution  equation 

-^u(L)  +  E(u(t),  A  )  =  0  (3.2) 

we  can  generate  the  basis  functions  using  solutions  to  the  stationary  problem  as  de¬ 
scribed  previously,  or  we  can  use  solutions  to  (3.2)  at  different  time  instances.  The 
reduced  order  solution  is  uR(t )  =  YL'iLi  ni(l )®i-  an(i  the  reduced  order  problem  is 


(—u  (I)  +  E{u  (/ ).  A),  <?.,) A-»xX  =  0,  i  =  1, 


,  m. 


We  performed  a  comparison  of  the  Lagrange  and  Hermite  basis  functions  for  the  sta¬ 
tionary  Navier- Stokes  equations.  Using  the  approach  described  above,  ie  computing 
solutions  at  various  Reynolds  numbers,  we  were  able  to  obtain  a  good  approximation 
using  3  to  5  basis  functions.  The  performance  of  the  Lagrange  versus  the  Hermite 
basis  functions  is  comparable.  Since  the  Lagrange  basis  functions  require  less  time  to 
compute  and  are  also  more  adaptable,  we  will  use  Lagrange  basis  functions  for  the  rest 
of  our  computations. 


3.1  Reduced  Order  Model  for  Viscous  Flows 

In  this  section  we  will  develop  a  reduced  order  model  for  the  unsteady  incompressible 
Navier-Stokes  equations.  The  reduced  basis  functions  will  be  generated  using  two 
different  techniques: 

1.  Solutions  to  the  time  dependent  Navier-Stokes  equations  at  various  time  in¬ 
stances, 

2.  Solutions  to  the  stationary  Navier-Stokes  equations  with  varying  boundary  ve¬ 
locity  control  force. 


6 


In  both  cases  we  will  use  Lagrange  interpolation  to  generate  the  reduced  basis  subspace. 
The  former  method  is  the  more  standard  approach  and  has  been  proven  effective  in  [10]. 
We  use  the  latter  approach  because  we  feel  it  leads  more  naturally  to  the  reduced  order 
control  problem.  Through  a  comparison  we  will  show  that  both  techniques  of  basis 
generation  produce  a  reduced  order  model  which  accurately  captures  the  dynamics  of 
the  time  dependent  Navier- Stokes  equations. 

We  will  first  illustrate  the  generation  of  the  reduced  basis  functions  from  the  time 
dependent  Navier  Slokes  equations.  Define  the  reduced  basis  subspace  as 

yR  =  span{ui  :  u;  =  u i  —  0, . . to} 

where  )  is  the  finite  element  approximation  at  time  T  to 

(4, vfe)  +  }-ea(*\vh)  +  b(u\  v*)  -  c(vh,ph)  =  (f,vh)  for  all  vA  G  Vq 

C(uh,qh)  =  0  for  all  qh  €  Sh 

(ufe,  v/l)r  =  (uft,  vfe)p  for  all  v^eVj 
(u(x,  0),  vA)  =  (u0(x),  vh)  for  all  vh  €  Vq. 

Define  the  reduced  order  solution  as 

m 

u'V)  =  4>0  +  »»  /(  /  )o,  , 

i= 1 

where  <po  =  u*  is  the  solution  at  the  final  time.  The  basis  functions  are 

defined  as 

4>i  =  U ■''(/,:)  -  Uh(li  |).  i  -  1,  .  .  TO. 

Thus,  cpo  has  boundary  conditions  as  the  inflow /out flow  conditions,  G  = 

V^  fl  Vq  have  homogeneous  boundary  conditions,  and  are  used  as  test  functions. 

We  will  now  illustrate  the  derivation  of  the  reduced  basis  functions  from  the  sta¬ 
tionary  Navier- Stokes  equations.  Define  the  reduced  basis  subspace  as 

VR  =  span{u{  =  :  u/l|rc  =  c,r,  i  =  1, . . . ,  to} 

where  the  finite  element  approximation  satisfies 

— a(uh,vh)  +  b(uh,uh,vh)  +  c(vhiPh)  =  (f,vft}  for  all  vA  G  Vq 
Re 

c(  uh,qh )  =  0  for  all  qh  G  Sh 

{  uh,  v/l)r  =  (u&,  vA)r  for  all  \h  G  Vq 

u|rc  =  I'iT 

where  Tc  is  the  control  portion  of  the  boundary,  Vi,  i  =  1. . .  .,m  are  the  boundary 
control  velocities,  and  r  is  a  unit  tangential  vector.  Define  the  reduced  order  solution 
as 

m 

U R(t)  =  <po  T  (/)<?, 

where  <po  =  uo  is  the  solution  with  zero  boundary  control  velocity,  and  the  basis 
functions  { <p., }''  ,  are  defined  as 

<pi  =  u m  —  a,-Uj  +  64-uq,  i  =  1 , . . . ,  to. 
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The  constants  a,  and  b,  are  chosen  to  force  homogeneous  boundary  conditions. 

Once  we  have  the  reduced  basis  subspace  V,'./  we  can  write  the  reduced  order 
equation  for  the  time  dependent  Navier-Stokes  equations  ( 2. 1 )-( 2.4)  as:  Find  u R(t)  G 
VR  such  that 

(uf ,  vfi)  +  ±-a{ uR,  v*)  +  b(uR ,  u*,  vfl)  =  (f ,  vR)  for  all  v*  €  V*  (3.3) 

(u(ar,  0),  vR)  =  (  u0 ( ■>' ) ,  v11)  for  all  vR  G  Vq  (  3.4) 

Note  that  the  divergence  free  term  does  not  appear  since  the  reduced  order  basis 
functions  satisfy  the  incompressibility  condition  (2.9).  Substituting  uR  =  <^0  +  YT=  1 4>i 
into  Equation  (3.3)  (3.1)  with  vR  =  <£j,  j  =  \. ...  .in.  we  obtain  the  following  system 
for  i  =  1, . . . ,  m: 

m  ^  171 


+  X  ^0  +  X  4')  =  (/:  &)  :  J  ~  1 . 

2  =  1  2  =  1 

771 

X  <A?)  =  (U0>  #X  J  =  1,  ■  ■  ■  ,  TO 

j=l 

We  can  write  this  system  as  the  following  system  of  nonlinear  ordinary  differential 
equations 

Ma(t )  +  A£(t)  +  col[Z(t)T VtZ(t)]™=1  -  /■'(/)  (3.5) 

Ma(0)  =  a0  (3.6) 

where  a  =  co/(m.  ...,  am)  €  D{  "  .  £  =  (1  a)T  G  Rm+1,  and  we  dehne  the  mass 
matrix  At  =  ( AT,y)  €  ]RmXm,  the  stiffness  matrix  A  =  (Aij)  G  JR,mX(m+1)j  the  forcing 
term  F  =  (Fj)  G  ]Rm,  and  the  initial  condition  U°  =  (E/f)  €  ]Rm  as  follows: 


t  j 

1— 1 

II 

II 

. . .  ,  TO, 

1—1 

i 

> 

1—1 

1 

> 

II 

i  =  2, . .  ,,m  +  1,  j  =  1, .  . 

•  ,  TO  +  1, 

■fi(0 

=  (/•'(/)•<•'/)•  '  -  I- 

•  •  •  ,  TO, 

=  (u0,O;),  i  -  1.  •  • 

. ,  m. 

The  matrices  Vt  G  ]R (m+1)x(m+1)  are  dehned  by 


(Vi)jk  =  (<f>j- 1  •  Vo/,  |.p,).  i  =  1, . . . ,  to,  j.  /,'  =  1, . . . ,  m  +  1. 

If  we  dehne  the  matrix  A"  as  ( ^  A"£ ),  =  ^  Pit,-  i  =  then  we  can  write 

Equation  ( 3.5)-( 3.6)  as 

Afd(i)  +  M£(2)  +  £{t)TN£(t}  =  F(t ) 


Afa(O)  =  U° 

This  is  a  dense,  low-order  system  of  nonlinear  ordinary  differential  equations  which  we 
solved  using  the  Implicit  Midpoint  rule 

A!'1'?  /  +  Atk+*  +  =  F*+s, 

Za  t  j  z 
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ak+1  =  2ak+^  -  ak. 

Ma°  =  U°. 

We  solved  this  quadratic  system  using  Newton’s  method.  Since  the  order  of  the  system 
is  small,  the  Newton  step  can  be  solved  efficiently  using  a  direct  solver. 

We  have  shown  that  a  reduced  order  model  can  be  developed  which  is  inexpensive 
to  solve.  In  the  next  section  we  will  show  that  this  model  can  accurately  capture  the 
original  dynamics  using  a  small  number  of  basis  functions. 

3.2  Computations  of  the  Reduced  Order  Model 

We  will  demonstrate  the  performance  of  the  reduced  order  model  for  the  time  depen¬ 
dent  Navier-Stokes  equations  on  the  forward  facing  channel  flow  problem.  As  described 
in  the  previous  section,  we  will  compare  two  techniques  for  the  selection  of  basis  func¬ 
tions.  The  basis  functions  generated  from  the  time  dependent  solution  are  shown  in 
Table  1. 


Basis  Function 

u0 

Ul 

u2 

u3 

u4 

Number  of  time  steps 

100 

200 

300 

400 

500 

Table  1:  Basis  Functions  Generated  from  the  Time  Dependent  Navier-Stokes  Equations 

The  basis  functions  are  defined  as  follows: 

0o  =  U4 

01  =  ui  -  u0 

0 2  =  U2  -  Ui 

03  =  U3  —  U2 

04  =  U4  -  U3. 

The  basis  functions  generated  from  the  stationary  solution  with  corresponding  bound¬ 
ary  velocity  are  shown  in  Table  2.  The  control  region  is  taken  to  be  the  portion  of  the 
lower  boundary  from  x  =  1  to  x  =  5. 


Basis  Function 

u0 

Ul 

u2 

u3 

u4 

U.5 

Boundary  velocity 

0 

0.05 

0.1 

0.15 

0.2 

0.25 

Table  2:  Basis  Functions  Generated  from  the  Stationary  Navier-Stokes  Equations 

The  basis  functions  are  defined  as 

0  o  =  u0 

01  =  u.5  -  f)U  I  +  4u0 

02  =  U.5  -  |u2  +  |u0 

03  =  U.5  -  |u3  +  |u0 

04  =  U.5  -  5U4  +  |u0. 
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In  Figure  2  we  see  the  full  solution  at  7  =  25.  The  solution  was  generated  from  the 
mixed  finite  element  approximation  described  in  Section  2.  In  Figure  3  we  see  the 
reduced  order  solution  using  basis  functions  generated  from  the  stationary  solution.  In 
Figure  4  we  see  the  reduced  order  solution  using  basis  functions  generated  from  the 
time  dependent  solution.  Since  we  used  the  solution  to  the  Stokes  equations  as  an 
initial  condition  for  the  full  solution,  we  will  use  the  projection  of  the  Stokes  solution 
onto  the  reduced  basis  space  as  an  initial  condition  for  the  reduced  order  model. 

Comparing  Figures  3  and  4  with  Figure  2,  we  see  that  the  reduced  order  models  do 
accurately  capture  the  dynamics  of  the  full  solution. 


Figure  2:  Full  Solution  at  T=25 


Figure  3:  Reduced  Order  Solution  with  Stationary  Basis  Functions  at  T=25 

For  a  closer  examination,  in  Figure  5  we  see  the  velocity  profiles  of  the  reduced  order 
models  and  the  full  solution  at  various  time  instances.  The  velocity  profiles  are  chosen 
at  x  =  1,2,  3, 4, 5.  These  were  chosen  because  this  is  where  the  recirculation  is  taking 
place,  and  it  is  this  property  which  we  want  to  capture  with  the  reduced  order  models. 
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Figure  4:  Reduced  Order  Solution  with  Time  Dependent  Basis  Functions  at  T=25 

4  Reduced  Order  Control  Problem 

In  this  section  we  will  develop  a  reduced  order  model  for  the  control  of  the  unsteady 
Navier- Stokes  equations.  From  this  reduced  order  system  we  will  compute  the  neces¬ 
sary  optimality  conditions  and  outline  a  method  of  solution  which  can  also  be  used 
to  compute  the  open  loop  control  for  control  problems  governed  by  partial  differential 
equations. 

We  begin  with  the  development  of  the  optimal  control  problem  for  viscous  incom¬ 
pressible  flows  with  boundary  velocity  control.  Let  us  formulate  the  optimal  control 
problem  as 


Minimize  77(u,  r(/Y)  =  j 

fT(H*)  + 
0 

‘^\v(t)\2)dt 

subject  to 

(w,  v'; 

)  +  a(u,  v)  +  b{  u,  u,  v)  -  c(  v,p) 

=  (f,v) 

for  all  v  6  Hq(R) 

h(u  ,q) 

=  0  for 

■  all  q  6  i2(fl) 

u 

=  uj  on  f 

u 

=  u(t)r 

on  rc 

(u(a?,  0),  v) 

o 

II 

),v)  for  all  v  €  Hg(fl) 

where  f  is  a  fixed  body  force,  m,  is  a  fixed  boundary  value  on  the  boundary  I  .  The 
function  u(t )  is  the  control  velocity  defined  on  the  control  portion  of  the  boundary  fc, 
and  r  is  a  unit  tangential  vector.  In  the  cost  functional,  /3  is  the  penalty  term  and  v(t ) 
is  the  control  force  with  the  relation  Jju(f)  =  v[t).  Note  that  by  including  v(t)  in  the 
cost  functional,  we  are  minimizing  the  rate  of  change  of  the  control  velocity  u(t).  The 
problem  we  will  be  studying  is  the  minimization  of  vorticity  in  some  subregion  of  the 
forward  facing  channel.  We  will  use 

^(uj  =  -^vorticity  =  /  J V  x  u|2  dx  (4.7) 

where  0*  C  R  is  the  subregion  where  we  wish  to  minimize  vorticity.  We  will  also  use 

-F(u)  —  -^tracking  —  lu  —  Udl  (4.8) 
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solid  line:  full  solution,  dashed  line:  reduced  order  solution  solid  line:  full  solution,  dashed  line:  reduced  order  solution 


Basis  functions  from  Stationary  Solution  Basis  functions  from  Time  Dependent  Solution 

Figure  5:  Velocity  Profiles  for  Forward  Facing  Channel  Flow 

where  is  the  solution  to  the  Stokes  equations.  Since  the  Stokes  equations  do  not 
have  a  recirculation  region,  we  are  minimizing  vorticity  by  tracking  to  this  solution. 

We  will  now  present  the  reduced  order  optimal  control  problem.  The  purpose  of 
designing  a  reduced  order  model  is  so  that  the  control  can  be  computed  in  an  efficient 
way,  thus  obtaining  a  suboptimal  control  law  for  the  full  system.  Define  the  reduced 
order  solution  as 

m  n 

XR(t)  =  (j)Q  +  X/'d/!''  + 

i= 1  j=l 

where  <f>o  corresponds  to  the  stationary  boundary  force,  ie  the  inflow/outflow  condi¬ 
tions,  {4>i}1£=1  are  test  functions  with  homogeneous  boundary  conditions,  and  1 

are  trial  functions  with  corresponding  boundary  control  force.  The  number  of  trial 
functions  n  will  be  determined  by  the  number  of  control  inputs  in  the  system.  The 
number  of  test  functions  m  should  be  chosen  large  enough  that  they  can  accurately 
capture  the  dynamics  of  the  system,  yet  small  enough  so  that  the  reduced  order  system 
does  not  become  ill-posed.  For  the  forward  facing  channel  flow  problem  with  tangen¬ 
tial  boundary  velocity  control,  we  have  one  trial  function  corresponding  to  the  control 
input,  and  we  obtained  good  results  using  2  to  4  test  functions. 
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We  can  write  the  reduced  order  controlled  unsteady  Navier- Stokes  equations  as 
Ma(t)  +  Ax(t)  +  x{t)Tj\fx(t )  +  bv{t )  =  0 
y*/(  / )  =  v{t). 

We  can  now  write  the  reduced  order  control  problem  as 

fT  i  a 

Minimize  J{x,v)  =  /  ( -xT  Qx  +  —  \v(t  )\2)dt 

Jo  2  2 

subject  to 

Ma(t)  +  Ax(t)  +  x{i)T  j\f  x(i)  +  bv(t)  =  0 
ii(t)  =  e(t), 

where  Qij  =  (V  X  (f>i,V  X  <f>j),  i,j  =  Note  that  we  can  write  this  control 

problem  in  a  general  form  as 

Minimize  »'(/))  -  /  (/(:<(/ ) )  f  li(  r(l)))<ll  (4-9) 

Jo 

subject  to 


4(0  =  /0(0)  +  (4.10) 

*(0|  =  so  (4.11) 

where  x(t)  =  |<i(  /)  u(t)]T  G  D{  "  is  the  state  and  v(t)  is  the  control.  If  we  let  h(v(l\)  = 
^|u(0|2,  then  we  can  write  the  necessary  conditions  of  optimality  as 

4(0  =  f(x(t))  +  Bv(t)  a.e.  in  (0, T),  s(0)  =  s0, 

-p( 0  =  fx(x(t))Tp(t)  +  lx(x(t))  a.e.  in  (0,T),  p(T)  =  0, 

v(t)  =  -^BTp(t)  a.e.  in  (0.7  ), 


which  is  the  well  known  Pontryagin  maximum  principle,  and  p(t)  G  i2(0,T)  is  the 
Lagrange  multiplier.  Thus  the  stale  coslate  pair  (x(t).p(t))  satisfies  the  two-point 
boundary  value  (TPBV)  problem 


4(1)  =  f(x(t))  -  i BBTp(t ) 

(4.12) 

p(  0  =  fx{x(t))Tp(t)  +  lx(x(t)) 

(4.13) 

with  s(0)  =  so  and  p(T)  =  0,  and  the  optimal  control  v(t)  for  ( 4.9)-( 4.11)  is  given  by 

00  =  -j-BTp(t). 

In  the  next  section  we  will  develop  a  solution  procedure  for  the  TPBV  problem  (4.12)- 
(4.13),  and  then  present  numerical  results  for  the  reduced  order  controlled  Navier- 
Stokes  equations. 
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4.1  Solution  Procedure  for  Reduced  Order  Control  Prob¬ 
lem 


In  this  section  we  will  describe  a  general  solution  procedure  to  solve  the  TPBV  problem 
(  1. 12)  (1. 13).  This  method  is  described  in  [10]  and  is  summarized  here  since  the  results 
are  used  directly  in  the  design  of  the  feedback  synthesis. 

We  will  discretize  (4.9)- (4.11)  using  the  Crari  k-Klcholson  method  for  the  dynamical 
constraint  (4.10),  and  the  trapezoidal  rule  for  the  cost  functional  (4.9),  instead  of 
directly  discretizing  the  TPBV  problem  (  1. 1 2)  (  1. 13).  That  is, 


N  . 

Minimize  JN  =  ^  (-(/(**  +  l(xk))  +  h{vk)) At 


subject  to 


k= 1 


xk  Af _1  =  +  Bvk,k  =  l,...,N 


(4.14) 


(4.15) 


where  N  At  =  T  and  x°  =  xq.  The  necessary  optimality  conditions  for  (4.14)-(4.15) 
are 


rpk  _  ^ 

At 

-Pk+1MPk  =  lUzk)T(Pk+Pk+1)  +  l^k) 


=  !(/(**}  +  /(**  '))  -  l~3BBTpk 


(4.16) 

(4.17) 


for  k  =  1  ,...,#  and  x°  =  a0,  pA+1  =  0,  and  the  optimal  control  to  (4.14)-(4.15)  is 
given  bv 

**  -  -i*V. 

As  a  consequence,  (4.16)-(4.17)  is  an  approximation  to  the  TPBV  problem  (4.12)- 
(4.13).  System  (4.16)-(4.17)  is  a  sparse  system  of  nonlinear  equations  in  xk ,  pk ,  k  = 
I , . . . ,  /V ,  which  can  be  solved  using  Newton’s  method: 


V+  --  !Jc  [B'(y,j]  /''(//,  ) 

where  y  =  co/(a1, . . . ,  xN  ^p1, . . .  ,pN ),  the  function  F(y)  =  /■'(./• 1 , . . . ,  aA  ,p1, . . .  ,pN ) 
is  given  by 


F(y)  = 


At 


-  |(/( X%  +  f(xk~1))  +  \BBTpk 


v-^-  +  \Uxk)T(pk  +  pk+1)  +  h{xk) 
and  the  Jacobian  has  the  following  matrix  structure 


F’  = 


A  C 
D  At 


where  A  is  block  lower  bi-diagonal,  and  C  and  D  are  block  diagonal  with  blocksize  m. 
The  diagonal  and  subdiagonal  blocks  of  A  are  given  by 

Ak,k  —  ~/\~f  2^x^X  ^  Ak+i,k  —  T  ^/a;(*£  )■ 
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The  diagonal  block  C  is  the  constant  block 

Ck,k  =  j3bbTi 

and  the  diagonal  block  D  is  given  by 

Dk,k  =  ^fxx{xk)T(pk  +  pk+1)  +  lxx(xk)- 

If  we  reorder  the  vector  y  =  co/(ad , ... ,  xN  ,px, . . .  ,pN)  to  y  =  col(xx ,  p1 ,  x2 ,  p2 , . .  . ,  xN ,pN ) 
then  the  Jacobian  is  a  block  tridiagonal  matrix  with  blocksize  2m.  Thus  the  Newton 
step  can  be  solved  using  a  block  tri diagonal  algorithm. 

4.2  Computations  of  the  Open  Loop  Control  Problem 

We  consider  the  control  problem  in  a  forward  facing  step  channel  flow.  The  goal  is  to 
minimize  the  recirculation  that  occurs  at  the  top  of  the  step  by  tangential  movement 
of  the  boundary  near  the  recirculation  region,  as  seen  in  the  schematic  in  Figure  6. 


U(t) 


Figure  6:  Schematic  for  Controlled  Forward  Facing  Channel  Flow 

The  basis  functions  {uj,  i  =  0, . . . ,  5}  were  chosen  as  solutions  to  the  stationary  Navier- 
Stokes  equations  with  boundary  control  velocity  as  shown  in  Table  3. 


Basis  Function 

u0 

Ul 

u2 

u3 

u4 

U5 

Boundary  Velocity 

0 

0.05 

0.1 

0.15 

0.2 

0.25 

Table  3:  Basis  Functions  Generated  by  Boundary  Velocities 
Define  the  reduced  order  solution  as 
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u  R(t)  = 

4>o  1  +  "CM 

i= i 

where  the  basis  functions  are  defined 

as 

4>o 

= 

u0 

4>i 

= 

u5  -  5ui  +  4u0 

<t>2 

= 

5  3 

Us  -  -u2  +  -u0 

<t>3 

= 

5  2 

us  -  -u3  +  -u0 

4>  4 

= 

5  1 

us  -  -U4  +  -Uo 

4>5 

= 

4(u5  -  u0) 

15 

Using  the  parameters  A t  =  0.05,  j3  —  ^ ,  and  an  initial  condition  r0  =  [0  0  0  0  0]r, 
we  obtained  convergence  of  the  two  point  boundary  value  problem  in  9  iterations.  The 
optimal  control  u(t)  is  shown  in  Figure  7.  The  flow  simulation  with  the  optimal  control 
input  is  shown  in  Figure  8. 


Figure  7:  Control  u(t )  for  Forward  Channel  Flow 


Figure  8:  Navier-Stokes  Equations  with  Control  Input 


5  Design  of  Feedback  Synthesis 

In  this  section  we  will  discuss  a  method  for  the  construction  of  a  feedback  synthesis 
v(t)  =  K(x(t))  based  on  the  reduced  order  control  system 

Ma[i)  +  Ax(t)  +  x{t)T  Af  x(t)  +  bv(t)  =  0 

u(t)  =  v(t) 

o(0  j  =  ctoi 
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(5.18) 

(5.19) 

(5.20) 


and  then  applied  to  the  unsteady  Navier- Stokes  equations. 

Consider  the  optimal  control  problem  (4.9)-(4.11)  in  IR  .  Then  the  optimal  control 
can  alternatively  be  given  by  the  feedback  form 

v(t.)  =  K(x{t.),t)  =  ^(x(t),Vx{x(t),t)), 

where  ).p(t ))  =  v(t)  minimizes  the  Hamiltonian 

P(0  ■  (/(*)  +  B(x)V)  +  h(v(t)). 

Here,  V  is  the  unique  viscosity  solution  to  the  Hamilton-Jaeobi-Bellman  (HJB)  equa¬ 
tion 

1  +  /(*)  ■  14  -  h*(-BT\ 4)  +  l(x)  =  0 

V(x,T)  =  0, 

where  B{x)  =  B  and  h*  is  the  conjugate  function  of  h.  If  we  use  h  =  j\v(t.)\2,  then  the 
conjugate  function  is  h*  =  —  \p{t)\2.  Thus,  the  optimal  feedback  solution  is  given  by 

v(  t)  =  ~^BTVx{x(  t),t)  =  A  (;/:( / ). /). 

Assuming  V(x(t),t)  G  C1  ( (0, 7  );  D{  "I,  we  have  the  following  relationship  between  the 
HJB  equation  and  the  Pont ry agin  maximum  principle 

P(t)  -  Vx(x 

where  the  pair  (x(t),p(t) |  satisfies  the  TPBV  problem  (4. 12)- (4. 13)  with  a-(0)  =  x0 
and  p(T)  =  0. 

Note  that  the  feedback  gain  I((x(t),t)  =  /\ '/  (:(  (/),  /  )  depends  on  the  final  time 
T.  Since  the  problem  under  consideration  is  autonomous,  we  can  write  KT(x(t),t )  = 
f\'l'  .(;/.•( /),  0).  Given  T  >  0,  define  the  function 

x0  —  $(x0)  =  Px°(  0), 

where  (x(t),p(t))Xo  is  the  solution  to  the  TPBV  problem  with  ar( 0 )  =  xq.  Then  we  can 
define  the  feedback  law 

Kt{x  o,0)  =  $(io,pIO(0))  =  -^Bt$(x  o) 

at  each  point  xq  G  !R  m.  Since  the  optimal  feedback  law  0)  is  impractical, 

we  replace  the  time  dependent  feedback  law  with  the  stationary  feedback  law 

v(  t)  =  K(x(t) |  =  KT(x(t),  0) 


for  T  >  0  sufficiently  large. 

Based  on  these  observations,  we  can  construct  a  suboptimal  feedback  law  by  car¬ 
rying  out  the  following  steps: 

1.  Choose  a  set  of  collocation  points  S  C  Era  and  calculate  the  solution  (x(t):p(t))Xo 
to  the  TPBV  problem  at  the  points  xq  G  S.  Thus  we  obtain  the  values  of  K  at 
the  points  iqGS. 
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2.  Use  an  interpolation  method  to  construct  the  feedback  synthesis  K ( x )  satisfying 

K(xq)  =  K(x o)  for  x0  G  S. 

For  Step  2,  we  used  an  interpolation  method  based  on  the  Green’s  function.  Assume 
we  have  the  set  of  collocation  points  {(a;o)i}f£i  G  £■  Define  the  feedback  law 

M 

K(x)  =  Y  G(x,  ( :C0),)  •  Vi 
i=  1 


with 

G((a0)i)  (•'•()  )j)  •  =  A'((a;0)«):  *  =  1,  ■  ■  ■ ,  M. 

If  the  Green’s  functions  G  satisfies  the  biharmonic  equation 

A2G(a;)  =  7(:r). 

then  G(a)  =  |.t|  for  n  =  3,  and  G(x)  =  x2(lnx  —  1)  for  n  =  2.  In  this  case,  ii'(a) 
minimizes 


j^n  |  Alt  (a)|2  dx  subject  to  K((x0)i)  =  K((x0){),  (xo)i  G  S- 

In  our  numerical  calculations,  we  used  a  family  of  Green’s  functions  of  the  form 

G(x,y)  =  \x  -  y\a,  1  <  a  <  4. 


5.1  Computational  Results  of  the  Feedback  Synthesis 

The  design  of  the  feedback  law  K(x)  requires  that  we  compute  the  solution  (x(t),p(t)j 
to  the  TPBV  problem  at  a  specified  number  of  collocation  points  (ao)«-  The  extent  of 
these  points  can  be  determined  from  the  solution  to  the  TPBV  problem  with  the  desired 
initial  condition.  In  Figure  9  we  show  the  reduced  order  solution  (aj(t),  yd/)) 

computed  from  the  TPBV  problem  using  the  Stokes  tracking  problem  (4.8)  with  an 
initial  condition  of  [0  0  0].  We  chose  the  collocation  points  uniformly  distributed  on 
the  intervals  [-2,2]  X  [—15, 15]  X  [—0.5,  0.5]  for  a  total  of  195  =  5  X  13  X  3  points  in 
the  ai,a25  u  directions,  respectively. 

We  used  the  Green’s  function  G(x,y )  =  \x  —  y\a  with  a  value  of  o  =  3.7.  In 
Figure  10  we  compare  the  interpolated  values  of  K(x)  against  the  exact  feedback  law 
K(x).  The  parameter  u  is  fixed  at  u  =  0.5,  and  the  values  of  K(x)  are  computed  at 
1281  uniformly  distributed  points  in  ( [—2,  2],  [—15, 15],  0.5)  €  If3.  We  see  that  this 
method  does  provide  an  accurate  interpolant.  In  our  calculations,  we  observed  that 
lower  values  of  o  provided  greater  fluctuations  in  the  interpolant. 

In  Figure  11  we  compare  the  values  of  the  open  loop  control  and  the  closed  loop 
control  for  the  reduced  order  control  system  (5.18)-( 5.20 ).  It  validates  the  feasibility 
of  our  approach. 

We  now  describe  how  to  apply  the  feedback  control  law  on  the  full  Navier-Stokes 
equations.  To  calculate  the  feedback  control,  we  project  the  full  solution  onto  the 
reduced  basis  space  by 

2 

Pu(/)  =  4>0  +  Y  “'('M  +  Mi'm, 


18 


Figure  9:  Reduced  Order  Parameters  from  Stokes  Tracking  Problem 


Figure  10:  Interpolation  of  BT<&(x)  with  u  fixed  at  0.5,  a  =  3.7 

where  (apt),  &%(t) )  is  determined  by  the  projection  equation 

2 

(“(*)>  fa)  =  Oo  +  'faWi  +  "( /)<>3 ,  j  =  1,  2.  (5.21) 

Equation  (5.21)  is  then  written  as 

Mot  =  co/[(u(i)  oo  u(t)(j)m+1: 4>j)],  j  =  1, 2.  (5.22) 

We  can  then  construct  the  feedback  law 

v(t)  =  /\'(<l  |  (/l.O.jf  /  l, '«(/)  ), 

where  («i(t),  02(h))  are  computed  from  Equation  (5.22).  In  Figure  12  we  show  the 
values  of  the  closed  loop  control  from  the  full  Navier-Stokes  equations  compared  with 
the  values  of  the  closed  loop  control  computed  from  the  reduced  order  model.  The 
open  loop  control  is  also  included  as  a  reference. 
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Figure  11:  Feedback  Control  versus  Open  Loop  Control 


Figure  12:  Feedback  Control 

5.2  Dynamic  Compensator  Design 

The  design  of  the  feedback  synthesis  requires  that  we  have  a  complete  observation 
of  the  state  at  each  time  step.  In  practice,  this  is  not  feasible.  For  the  full  Navier- 
Stokes  equations,  rather  than  projecting  the  solution  onto  the  reduced  basis  space 
and  computing  the  feedback,  we  would  like  to  be  able  to  compute  the  feedback  input 
v(t.)  using  a  partial  observation  of  the  state.  In  this  section,  we  propose  a  method  of 
constructing  a  state  estimator  to  complete  the  compensator  dynamics. 

The  reduced  order  controlled  Navier- Stokes  equations  can  be  written  as 

dx  ^ 

-  +  £(*,„)  =  0 

z(0)  =  x0, 

where  v  =  K(x)  is  the  feedback  control  law  and  x  =  [a  u]T  are  the  parameters  of  the 
state.  Assume  we  have  a  partial  observation  of  the  full  Navier-Stokes  equations 

y(t)  =  h(u(t)), 
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where  h  is  a  vector- valued,  bounded,  linear  functional.  We  model  the  state  estimator 
as 

%  +  E(x,v)  =  G(y  -  Hit) 
at 

®(0)  =  *o, 

where  v  =  Ii(x)  is  the  feedback  control  law,  x  is  the  state  estimator,  and  Hj  = 
h(<i>j).  j  =  \ ...  ../m  is  the  matrix  representation  of  the  observation  functional  h  on  NR. 
We  compute  the  observer  gain  G  using  the  extended  Kalman  filter.  Let  A  =  F(x.  r) 
be  the  linearized  map  at  x.  Then  the  gain  G  is  given  by 

G  =  EMTiT15 

where  the  error  covariance  E  is  the  solution  to  the  Ricatti  equation 

4-E  =  AS  +  HA  -  YiHt R-1  HTi  +  Q, 

at 

and  Q ,  R  are  the  covariance  matrices  of  the  system  and  observation  noise,  respectively. 


6  Conclusions 

We  have  presented  a  reduced  order  model  based  on  the  reduced  basis  method  that  can 
be  used  for  the  simulation  and  control  of  viscous  incompressible  flow  problems.  This 
reduced  order  model  is  inexpensive  to  solve  and  accurately  captures  the  dynamics  of  the 
original  system  using  a  small  number  of  basis  functions.  The  feasibility  of  the  reduced 
order  model  for  flow  control  was  demonstrated  on  a  boundary  velocity  control  problem. 
A  low  order  feedback  control  law  was  developed  and  successfully  implemented  on  the 
full  Navier-Stokes  equations.  Through  our  numerical  experiments,  the  feasibility  of 
this  feedback  law  was  demonstrated  for  the  control  of  viscous  flows.  With  the  state 
estimator,  we  have  a  complete  compensator  design  for  the  full  Navier-Stokes  equations. 
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